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Introduction 


Hydrodynamic  predictions  of  waves  can  be  made  via  numerical  solution  of  the  Navier-Stokes  equations 
subject  to  the  appropriate  free-surface  boundary  conditions.  For  practical  flow  of  naval  hydrodynamic 
interest,  simplifications  to  these  equations  must  be  made  in  order  to  make  the  computations  feasible.  Two 
approaches  for  this  exist:  potential  flow,  where  viscous  and  rotational  effects  are  formally  eliminated 
from  the  equations,  and  Reynolds -averaged  approaches,  where  the  equations  are  time-averaged.  Potential 
flow  calculations  are  significantly  faster,  but  cannot  account  for  viscous  effects  or  turbulence.  Reynolds- 
averaged  Navier-Stokes  (RANS)  calculations  require  significantly  more  computational  effort,  as  well  as 
turbulence  modeling,  but  constitute  a  framework  which  can  accommodate  both  viscous  and  turbulence 
effects.  The  research  described  herein  is  motivated  by  the  need  for  wave  breaking  models  in  Reynolds- 
averaged  Navier-Stokes  (RANS)  calculations. 

Although  recent  developments  in  level-set  and  volume-of-fluid  numerical  methods  seem  to  have  relaxed 
the  need  for  detailed  wave-breaking  models  in  the  RANS  predictions  in  naval  hydrodynamics,  the 
absence  of  RANS  model  for  breaking  waves  means  some  physical  processes  are  not  captured.  The 
generation  of  free-surface  roughness  and  surface  currents,  the  net  dissipation  of  wave  energy  during 
breaking  or  the  production  of  bubbles  are  examples  of  processes  not  captured  by  standard  RANS 
modeling.  Numerical  stability  (via  numerical  dissipation)  and  geometric  generality  (via  surface  capturing) 
alone  cannot  fully  replace  the  need  for  models  of  the  unresolved  scales,  though  they  will  allow  for  some 
level  of  successful  calculations. 

An  approach  for  developing  RANS  models  for  breaking  waves  can  be  outlined  as  follows.  Given  a 
particular  numerical  methodology  in  a  predictive  engineering  application  (VoF,  Level  set,  potential  flow, 
ect.),  the  choice  of  which  fixes  what  can  be  represented  in  the  computation,  a  complementary  averaging 
(or  filtering)  process  must  defined  and  a  new  set  of  ‘RANS’  equations  developed.  While  the  resultant 
equations  may  appear  the  same,  the  models  needed  may  not  be.  The  unresolved  fluid  dynamics  in  a  VoF 
calculation  are  fundamentally  different  from  the  unresolved  fluid  dynamics  in  a  level-set  calculation. 
Given  the  ‘RANS’  equations,  the  filtering  process  and  the  relevant  experimental  and  numerical  data, 
develop  the  functional  form  for  the  model  of  the  unknown  terms  relevant  to  breaking  and  surface 
roughness.  The  work  herein  was  intended  to  provide,  when  fully  executed,  data  for  use  in  the  model 
development. 


Research  objectives 

The  primary  goal  of  this  research  effort  was  the  development  and  application  of  a  numerical  capability 
for  the  simulation  of  breaking  waves  including  the  effects  of  fluid  viscosity,  surface  tension  and  the 
generation  of  surface  roughness  and  turbulence.  The  capability  sought  was  to  maximize  fidelity  to  the 
physics  of  the  free-surface,  requiring  the  numerics  to  predict  all  relevant  scales. 

The  intent  was  to  develop  the  ability  to  both  simulate  the  breaking  event  and  to  have  a  platform  in 
which  models  could  be  developed  and  tested.  A  fully  coupled  formulation  was  adopted  in  the 
development  of  the  numerical  algorithm.  The  primary  benefits  of  this  approach  are  the  satisfaction  of 
mass  conservation  and  the  numerically  implicit  enforcement  of  strongly  coupled  boundary  conditions. 
The  algorithm  developed  is  remarkable  robust,  while  maintaining  high  order  of  accuracy.  The  algorithm 
was  implemented  in  two  and  three  dimensions  and  several  breaking  wave  simulations  were  performed. 
The  algorithm  was  developed  to  allow  both  model-free  DNS-like  calculations  and  RANS  calculations. 
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An  additional  objective  was  to  develop  an  solution  scheme  which  allowed  the  free-surface  to  become 
multiply  defined.  The  intent  of  this  adaptation  was  not  to  solve  for  plunging  wave,  but  rather  to  simulate 
the  roughness  found  in  less  aggressive  breaking  waves.  An  appropriate  solution  to  this  objective  was 
never  achieved. 

The  algorithm  was  applied  to  a  variety  of  conditions,  with  mixed  results.  The  evolution  of  small  scale 
two  dimensional  waves  could  be  followed  fully  through  breaking  and  into  the  decay  state.  Longer,  two- 
dimensional  waves  have  also  been  simulated,  but  without  sufficiently  high  resolution,  the  breaking 
process  is  insufficiently  dissipative.  Many  of  these  two-dimensional  simulations  become  unphysical  and 
fail.  Three-dimensional  waves  were  also  simulated,  but  the  algorithm  was  insufficiently  vetted  to  perform 
the  desired  full  three-dimensional  breaking  waves. 

Approach 

The  numerical  approach  adopted  for  the  study  of  breaking  waves  was  to  develop  a  robust  and  general 
solution  technique  for  the  Navier-stokes  equations  including  the  full  non-linear  boundary  conditions. 

Since  the  objective  was  the  development  of  numerical  dataset  of  simulated  breaking  waves  for  use  in  the 
development  of  engineering  models,  the  numerical  approach  needed  to  be  able  fully  represent  the  physics 
of  the  breaking  waves. 

Several  approaches  were  excluded  from  the  list  of  potential  numerical  methods.  The  excluded  numerical 
methods  included  Level-Set  methods,  Volume-of-Fluid  methods  and  any  approach  requiring  the 
linearization  of  the  free  surface.  While  level-set  methods  can  be  constructed  to  be  very  robust,  in  general 
the  only  satisfy  the  free-surface  boundary  conditions  in  a  weak  sense.  The  free  surface  boundary  is 
represented  as  a  distribution,  with  compact  support,  in  the  vicinity  of  the  mean  free  surface.  As  a  result, 
some  physical  processes  may  not  be  fully  represented.  On  the  other  hand,  VoF  methods  may  employ  a 
discrete  representation  of  the  fluid  volume  and  a  second  representation  for  the  interface,  such  that  the 
union  of  these  grids  represents  the  free  surface  and  the  fluid  volume.  However,  the  length  scales 
(boundary  layer  thickness)  normal  to  the  air- water  interface  scales  as 

I  ~  Re“1/2,  (1) 

and  fully  resolving  this  interface  is  not  practical  using  the  uniformly  spaced  grid  often  employed  in  VoF 
methods.  This  does  not  imply  that  these  methods  are  not  useful  in  many  free-surface  problems,  but  rather 
they  may  not  capture  the  physical  processes  needed  in  the  development  of  RANS  like  models  for  the 
interface. 

The  approach  adopted  was  a  surface  following  time  dependent  structured  grid,  a  high  order  compact  finite 
difference  representation  and  a  Krylov  space  algorithm  for  solving  the  fully  coupled  Navier-Stokes 
equations.  This  representation  satisfies  the  conservation  of  mass  locally  and  globally  on  a  fully  deforming 
grid.  The  primary  limitation  of  this  approach,  aside  from  computationally  efficiency,  is  the  inability  to 
generate  usable  grids  when  the  surface  slope  exceed  vertical. 
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Governing  Equations 

Navier-Stokes  (RANS)  governing  equations  are: 

V  ■  u  =  0  (2) 

*  +  2 ((u  V)u+  V .  (uu))=  -V(p+  gy)+  V  ■  ((v  +  v,)Vu)+  Fb  (3) 

where  Fb  is  an  body  force  used  to  model  the  forcing  of  submerged  foils  and  Ve  is  an  eddy  viscosity, 
employed  as  needed.  Note  the  hydrostatic  pressure  is  retained  in  the  governing  equations,  rather  than 
being  included  in  the  boundary  condition.  Boundary  conditions  imposed  on  the  free  surface  are  the 
kinematic  boundary  condition 


— -  +  u  V77  =  v  on  y  =  ;Xx,z,t),  (4) 

at 

and  the  dynamic  boundary  conditions 

n  a  nT  =-Pa-^V  n 

tx  C  nT  =  tx  xa  (5) 

tz  c  nT  =  tz  Ta 

where  Pa  and  Ta  are  the  normal  and  shear  stress  imposed  by  the  atmosphere.  na,  tx  and  tz  are  the  free 
surface  normal  and  tangential  vectors  ,  respectively,  and  O  is  the  stress  tensor: 

o  =  -pI  +  i(Vu+(Vu)T). 

Either  shear-free  and  no-slip  boundary  conditions  can  be  employed  on  the  bottom  of  the  domain  which 
may  be  an  arbitrary  shape. 


Numerical  Approach 

The  time  stepping  algorithm  is  a  mild  modification  to  the  common  Adams-Bashforth/Crank-Nicholson 
scheme.  In  such  a  scheme,  the  convective  terms  as  advanced  in  time  explicitly  by  an  Adams-Bashforth 
representation.  The  linear,  but  coupled  viscous  and  pressure  terms  in  the  NSE  are  advanced  implicitly 
with  the  Crank-Nicholson  scheme.  While  the  Crank-Nicholson  scheme  is  unconditionally  stable  for 
viscous  flows,  as  the  Reynolds  number  increases,  the  pressure  field  will  develop  a  ‘checker  board 
instability’  in  time.  Therefore  to  prevent  the  formation  of  the  temporal  instability,  the  Crank-Nicholson 
scheme  is  replace  by  the  so-called  theta  scheme: 

df  „gFn  +  (l-^)fn_1 
dt  n  At 
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where  0  <  0  <  1.  When  9  =  0.5,  the  scheme  is  standard  Crank-Nicholson,  while  in  this  work  a  value  of 
0  =  0.51,  was  usually  sufficient  to  prevent  oscillations  in  time.  For  simplicity,  the  description  of  the  time 
advancement  scheme  assumes  the  Crank-Nicholson  scheme 


The  outline  of  the  time  advancements  algorithm  is  as  follows: 

„■  +  *  [V p]"  [V  ■  ((v  +  v,)Vu)J  =  u-1  + 1  [Vpf  -  T  [v  ■  ((r  +  v,)Vu)' 

+  a[(|N-1-iN-J)+(Fb-ge2); 


V  ■  u11  =  0 


„  dr\ 

u  — - 

-^Vn  =77n-1  +  4V“1-— u 

drj 

\dx 

2  2  2 

_dx 

where 

N  =  i((u  V)u+V  (uu)) 


(6) 


is  the  nonlinear  term,  and  u  is  an  extrapolated  value  of  the  surface  velocity  at  forward  half  step  in  time. 
The  computations  are  performed  in  the  physical  domain  or  grid  (the  velocities  remain  the  usual  Cartesian 
velocities),  but  the  derivative  are  all  evaluated  on  an  orthogonal  computational  grid.  The  mapping 
between  the  physical  grid  and  computational  grid,  represented  by  the  metrics  is  general  and  grid  motion  is 
fully  accounted  for. 


Collectively  these  equations  are  a  set  of  3  2-D  plus  1  1-D  coupled  equations  in  two  dimensions  or  4  3-D 
plus  1  2-D  coupled  equations  in  three  dimensions.  These  equations  can  be  written  formally  as 


An+1 

Gn+1 

B"+1 

un+1 

rn 

bn 

u 

>i 

u 

u 

Dn+1 

0 

0 

Pn+1 

= 

0 

+ 

0 

En+1 

0 

An+1 

J 

T]n+l 

1 

1 _ 

K 

(7) 


The  first  row  represents  the  momentum  equation  and  the  second  is  the  continuity  equation.  The  third  row 
is  the  free-surface  kinematic  boundary  condition.  Most  of  the  entries  in  the  left  hand  side  are  self-evident, 
but  two  are  worthy  of  note.  In  the  top  row,  B  ”+1  represents  the  coupling  of  the  free-surface  elevation  into 

the  free-surface  dynamic  boundary  conditions,  while  in  the  third  row,  En+1,  is  the  coupling  the  surface 
velocity  into  the  kinematic  free  surface  boundary  condition.  The  metrics  are  included  in  the  terms  of  the 
matrix,  which  must  be  updated  each  time  step  since  the  grid  metrics  are  time  dependent. 

The  common  approach  in  solving  the  incompressible  Navier-Stokes  equations  is  to  decompose  or  split 
them  into  two  or  more  sequential  processes  or  operators.  The  operator  splitting  is  done  to  decouple  the 
viscous  and  convective  processes  from  the  calculation  of  the  pressure  and  enforcement  of  continuity. 
While  robust  algorithms  can  be  developed  by  employing  splitting,  when  operator  split  algorithms  are  used 
in  Direct  Numerical  Simulations,  the  continuity  constraint  may  not  be  completely  satisfied  at  the 
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boundaries.  Furthermore,  the  coupling  exhibited  in  the  free-surface  boundary  conditions  cannot  be 
handled  consistently  in  a  split  algorithm. 

The  alternate  and  adopted  approach  is  to  retain  the  fully  coupled  system  of  equations,  and  implementing 
an  algorithm  which  allows  the  satisfaction  of  the  boundary  conditions  and  the  continuity  constraint  in  a 
numerically  semi-implicit  form  time  stepping  algorithm.  This  coupled  set  of  algebraic  equations  must  be 
solved  at  each  time  step.  To  keep  the  problem  manageable,  an  algorithm,  which  does  not  require  the 
explicit  construction  of  the  algebraic  system  is  necessary. 

A  matrix  free  pre-conditioned  bi-conjugate  gradient  scheme  (BiCGstab)  was  employed  to  solve  the 
system  of  equations.  The  BiCGstab  schemes  have  been  shown  effective  in  algebraic  systems  which  are 
strongly  asymmetric  and  not  positive  definite,  as  long  as  an  appropriate  pre-conditioner  can  be  developed. 
To  be  effective  in  this  calculation,  the  pre-conditioner,  which  is  a  computationally  inexpensive 
approximation  to  the  system  above,  must  act  as  a  regularizing  process.  An  incomplete  Uzawa-scheme  was 
used  as  the  preconditioner.  The  Uzawa  scheme  is  essentially  an  approximate  LU  matrix  decomposition  of 
the  equation  (7).  The  incomplete  LU  decomposition  as  applied  to  the  equation  above  is  similar  to  the 
operator  splitting  scheme  noted  above.  The  preconditioner  is  applied,  in  parallel,  on  the  decomposed 
domain.  The  BiCGstab  algorithm  provides  the  global  coupling. 

For  sufficiently  smooth  grid,  such  as  that  in  Figure  1,  the  iteration  scheme  forces  the  residual  of  equation 
(7)  to  ~10  13  in  150  to  250  iterations.  Even  for  very  ill-constructed  unsmooth  grid,  such  as  those  observed 
when  the  simulated  breaking  event  becomes  unphysical,  the  numerical  scheme  will  continue  to  converge, 
albiet  slowly,  to  a  (unphysical)  solution.  While  computing  unphysical  solution  is  the  objective,  a  robust 
solution  technique  was  required. 


Geometry  and  Grid  Generation 

The  formulation  of  the  grid-transformed  equations,  described  above,  is  completely  general.  The 
computational  grid  is  not  required  to  be  orthogonal  within  the  domain  or  at  the  free  surface.  The  grid  may 
include  a  multiply  defined  free-surface.  The  only  limitation  as  formulated  was  that  the  domain  by  simply 
connected:  re-entrant  breaking  waves  wherein  enclosed  air  cavities  occur  are  excluded.  Although  surface 
following  algorithms  may  best  predicted  capillarity  dominated  wave  physics,  are  ill-suited  for  re-entrant 
wave  breaking.  An  objective  of  the  research  was  to  develop  a  grid  generation  scheme  allowing  for 
multiply  defined  surfaces.  This  objective  was  to  address  the  occurrence  of  highly  contorted  surfaces 
where  in  the  slope  may  exceed  vertical,  but  not  include  reentrant  flows,  but  a  satisfactory  method  of 
generate  a  smooth  grid  for  the  multiply  defined  geometry  was  never  developed 

The  waves  to  be  simulated  were  assumed  from  the  start  to  exist  in  a  periodic  domain  so  that  the 
emphasis  could  be  placed  on  the  development  of  the  numerical  method.  Two  approaches  for  grid 
generation  were  considered:  algebraic  grid  generation  and  hyperbolic  grid  generation.  For  sufficiently 
benign  free-surfaces  the  hyperbolic  grids  were  acceptable,  though  more  costly  than  simple  algebraic 
grids.  However  generating  a  smooth  hyperbolic  grid  when  the  free-surface  becomes  sufficiently  distorted 
is  difficult.  The  easier,  though  more  restrictive,  approach  was  to  use  algebraic  grid  generation.  The 
algebraic  grids  are  intrinsically  smoother  and  resulted  in  better  convergence.  The  algebraic  grids  were 
generated  via  a  spline  technique,  which  allows  for  easy  adjustment  to  keep  the  grid  nearly  orthogonal  the 
free  surface.  Figure  1  displays  an  example  of  the  streamwise  velocity  of  an  artificially  steep 
wave  (kct  =  1,  X  =  1  m) ,  along  with  the  grid. 


7 


Enhanced  Numerical  Modeling  of  Breaking  Waves  (2). doc 


Figure  1  Sample  wave  calculation  showing  streamwise  velocity  and  deforming  grid. 


Figure  2  displays  the  streamwise  velocity  over  a  shear  free  hump,  along  with  an  instance  of  the  grid. 
The  (periodic)  channel  is  a  meter  long  and  10  cm  deep.  The  objective  was  to  examine  and  model  the 
production  of  turbulence  in  a  breaking  zone,  modeled  as  a  hydraulic  jump.  A  multiblock  boundary  exists 
at  y  =  0.5.  As  the  computation  proceeded,  the  face  of  the  hydraulic  jump  approach  vertical  and  the 
calculation  failed.  Mass  was  conserved  during  the  calculation. 


0.05 
-0.00 
-0-05 

-0.10 

-0.15 

Figure  2.  Streamwise  velocity,  flowing  to  the  right,  over  a  shear-free  hump.  Mean  velocity  is  0.5  m/s. 

Slow  velocity  is  purple,  faster  velocities  are  red. 

As  implied  above,  the  domain  was  decomposed  in  several  block  and  a  parallel  algorithm  was 
implemented  using  the  MPI  standard. 

Extensions 

Two  extensions  to  the  two  dimensional  algorithm  were  to  be  implemented.  As  the  calculations 
proceeded,  turbulence  models  were  to  be  implemented.  These  turbulence  models  were  to  enhance  the 
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dissipation  of  wave  energy  and  control  the  breaking  process.  The  second  extension  was  to  implement  the 
algorithm  in  three  dimensions. 

Turbulence  Modeling 

Two  types  of  turbulence  models  were  implemented:  a  simple  Smagorinsky  type  LES  model  and  the 
standard  k-e  two  equation  model.  The  former  was  a  simple  add-on  sometimes  used  in  initializing  the 
code.  The  two-equation  model  was  exercised,  but  never  thoroughly  vetted.  The  primary  missing  element 
in  the  model  was  a  term  to  initiate  the  production  of  turbulence  kinetic  energy  when  the  spilling  starts. 
Several  were  considered  based  on  the  local  surface  compression  on  the  face  of  the  wave,  but  none  were 
successfully  developed  in  the  time  allocated. 

Extension  to  Three  Dimensions 

The  algorithm  was  also  implemented  in  three  dimensions.  The  3-D  code  exhibits  the  same  iterative 
convergence  characteristics  as  its  two-dimensional  counterpart.  The  code  retains  all  the  same  grid  issues 
and  numerical  approaches  as  the  two-dimensional  version.  Some  limited  results  will  be  presented  below 
(See  Figure  6.) 

The  code  was  not  sufficiently  vetted  before  the  end  of  the  contract  to  effectively  use  the  HPC  resources 
available  at  the  time  and  is  too  computationally  intensive  for  the  currently  available  resources.  Tuning  the 
inter-process  communication  was  accomplished  for  the  two  dimensional  application,  but  remained 
inefficient  in  the  three-dimensional  code. 


Two  Dimensional  Calculations 

The  role  of  the  two-dimensional  calculations  was  for  validation  and  optimization  of  the  algorithm  in  the 
limit  of  strong  small  scale  dynamics.  While  these  simulations  may  be  of  questionable  utility  in  the  context 
of  naval  hydrodynamics,  they  are  a  necessary  step  towards  developing  tools  for  predicting  the  larger  scale 
breaking. 


‘Wind  Driven’  Breaking 

A  set  of  small  scale  breaking  wave  simulations,  referred  to  as  wind  driven,  were  performed  for  the 
purpose  of  (1)  understanding  the  limits  and  limitations  of  the  algorithm,  (2)  testing  and  optimizing  the 
computer  code  and  (3)  to  provide  insight  into  the  process  by  which  turbulence  and  surface  roughness  are 
generated.  These  simulations  calculate  the  evolution  of  an  initially  linear  wave  driven  to  breaking  by  a 
phase-locked  external  ‘atmospheric’  pressure.  Such  a  calculation  is  presented  in  Figures  3-5,  wherein  a 
small-scale  wave  is  driven  past  breaking.  In  Figure  3,  the  vorticity  is  generated  at  the  interface  near 
x=0.625.  The  source  of  the  vorticity  is  driven  downstream  by  the  wave  propagation  and  as  expected  a 
packet  of  resonant  capillary  waves  are  observed  ahead  of  this  point.  (The  decaying  vortex  seen  upstream 
is  a  consequence  of  the  periodic  domain.)  The  vorticity  is  unsteady,  rolling  up  into  nearly  discrete 
vortices.  However,  the  generation  of  vorticity  occurs  predominately  at  the  localized  source. 
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Figure  3  Tangential  vorticity  in  late  breaking  of  a  10  cm  wave.  The  vorticity  ranges  from  -3250  sec-1  to 
3750  sec-1.  The  frequency  of  the  wave  is  approximately  25  rad  sec-1.  The  vertical  bar  at  x=0.5  indicates 
the  block  boundary. 

The  surface  elevation  characteristics  of  this  wave  are  presented  in  Figure  4.  The  salient  features  are  (1) 
the  strong  depression,  and  resulting  curvature  near  x=0.35  and  (2)  the  strong,  but  compact,  negative  peak 
in  curvature  at  the  source  vorticity.  The  magnitude  of  the  instantaneous  slopes  can  exceed  unity. 

The  complicated  interaction  of  the  compact  vorticity  and  the  strong  depressions  like  that  occurring  at 
x=0.35  may  be  a  process  by  which  air  is  entrained  at  the  interface.  When  simulations  of  this  type  are 
performed  at  lower  resolution,  these  interaction  are  not  as  intense.  With  increasing  resolution,  these 
cavities  have  been  observed  to  grow.  Due  to  grid  limitations,  these  cavities  cannot  pinch  off  and  the 
growth  leads  to  unphysical  results.  To  prevent  the  uncontrolled  growth  of  these  cavities,  a  special  sixth- 
order  filter  is  applied  to  the  surface  elevation. 


Figure  4  Surface  elevation,  slope  and  curvature  of  wave  in  Figure  3. 
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The  plot  of  the  curvature  in  Figure  4,  exhibits  small  wiggles,  which  are  probably  erroneous.  Once  the 
wiggles  are  created,  real  or  otherwise,  they  propagate  like  capillary  waves,  with  a  phase  speed 
proportional  to  the  reciprocal  of  the  square  root  of  their  wavelength.  Simulating  these  small  capillary 
waves  requires  an  implicit  numerical  scheme  to  reduce  the  time  stepping  constraint  in  explicit  schemes. 

A  final  observation  relevant  to  the  development  of  wave  breaking  models  can  be  made.  In  Figure  5,  the 
fluid  velocity  at  the  surface  is  resolved  into  the  surface  normal  and  tangential  directions.  Note  the  strong 
negative  rate-of-strain  where  the  vorticity  is  generated,  while  the  rate-of-strain  has  a  dipole  character  near 
the  depression.  As  a  result  the  will  be  a  flux  of  vorticity  into  the  fluid  in  the  former  case,  but  not  in  the 
later  case.  This  observation  is  important  in  the  context  of  the  development  of  RANS  models  for  breaking 
waves.  In  addition  to  providing  some  unsteadiness,  this  region  for  the  production  or  flux  of  vorticity 
creates  the  shear-layer  riding  on  the  face  of  the  breaking  wave. 


Figure  5.  Fluid  velocity  at  surface  of  the  wave  in  Figure  3  resolved  into  surface  normal  and  tangential 
components.  .  For  comparison,  the  phase  velocity  of  the  linear  wave  in  these  non-dimensional  units  is 
approximately  0.4.  The  surface  rate-of-strain  in  the  bottom  frame. 
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Three-Dimensional  Calculations 


A  preliminary  and  limited  set  of  three-dimensional  wave  simulations  for  the  purpose  of  testing  and 
validating  the  code  were  performed.  Two  geometries  were  considered:  interacting  oblique  waves  and 
recirculating  open  channel  flow  with  a  bottom  bump.  These  cases  validated  the  general  extension  of  the 
numerics  to  three  dimensions.  Additionally,  the  inter-block  communications  used  in  the  domain 
decomposition  were  tested  and  validated. 

A  sample  calculation  of  recirculating  flow  over  a  shear  free  hump  is  shown  in  Figure  6  &  7.  This 
computational  domain  consists  of  six  blocks.  The  three  frames  showing  unsteady  pressure  in  Figure  6  are 
separated  in  time  by  the  H/U  where  H  is  the  channel  depth  and  U  is  the  mean  velocity.  Figure  7  is  the 
streamwise  velocity  at  the  last  frame.  The  pressure  field  drops  rapidly  behind  the  bump  resulting  in  a 
strong  surface  depression.  As  a  result  of  Bernoulli’s  law,  the  velocity  above  the  hump  is  also  increased. 
As  the  calculation  proceeds,  the  flow  attempts  to  develop  a  hydraulic  jump,  but  the  simulation  was  never 
performed  with  sufficient  resolution. 


Figure  6.  The  unsteady  pressure  in  a  recirculating  channel.  The  frames  are  separated  by  H/U  in  time. 
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Figure  7..  Streamwise  flow  over  shear-free  hump.  Flow  is  to  the  right. 


Summary 

The  primary  objectives  were  the  development  and  application  of  the  capability  to  simulated  breaking 
waves.  This  was  seen  as  part  of  the  larger  and  longer  term  objective  of  developing  RANS  models  for  the 
effects  breaking  waves. 

Algorithms  were  developed  for  simulating  waves  in  two  and  three  dimensions.  To  a  limited  extent  these 
computer  codes  are  capable  of  meeting  the  objective  of  simulating  breaking  waves.  The  numerical 
algorithms  are  robust  in  both  two  and  three  dimension,  but  critical  limitations  remain. 

In  two-dimensions,  small  scale  unsteady  breaking  was  simulated.  The  development  of  surface 
roughness  and  the  possible  entrainment  of  air  were  observed.  In  some  strongly  forced  open  channel 
simulations,  the  surface  could  become  sufficiently  energetic  that  possible  fluid  ejection  was  observed. 
With  increase  wavelength,  the  restraining  influence  of  surface  tension  was  reduced  and  the  simulations 
become  more  difficult.  This  is  more  the  result  of  the  limited  grid  implementation,  than  in  the  algorithm 
itself. 

A  similar  algorithm  was  implemented  in  three-dimensions,  but  insufficient  time  was  allocated  to  its 
development  and  validation.  While  the  capabilities  are  the  same  as  its  two-dimensional  counter  part, 
substantially  greater  resources  are  required. 
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